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I. VELOCITY AUTOCORRELATION II. KINETIC ENERGY DISTRIBUTIONS 

FUNCTIONS 

Figure [2] presents the distribution of the kinetic ener- 
Figure [T] presents the transverse velocity autocorrela- gies among the particles in the system, 
tion functions for different damping strengths, in the qui¬ 
escent system. 






FIG. 1. Transverse velocity autocorrelation functions for 
particle mass m = 1 and different damping magnitudes C, 
in (top) quiescent systems, i.e., Ei = 0, at To = 0.16 and 
(bottom) sheared systems, with Ei and To as indicated in 
the legend. Note that the values of To used for the flows at 
C = 1 were chosen to match the kinetic temperatures of the 
strongly underdamped flows (^ = 10“®) at the same rescaled 
shear rate Ei. The dashed black line is an envelope fit with 
the function exp) ^ ^ ) 
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FIG. 2. Distribution of the kinetic energies ck = ^rnv^ of 
the individual particles in the athermal and thermostatted 
systems (as indicated in the legend). Solid black lines are fits 
to Boltzmann distributions oc 


III. ROBUSTNESS OF OUR FINDINGS: TEST OF 
A DIFFERENT THERMOSTAT 


To test the robustness of our results, we subsitute 
a Langevin thermostat for the Dissipative Particle Dy¬ 
namics scheme used in the main text. Figure [3] proves 
that, for low enough damping strength, the athermal 
Langevin flow curve also becomes nonmonotonic, with 
a rate-weakening region. More importantly, the ther¬ 
mostatting procedure described in the main text is also 
operational here, as evidenced in the figure. This con¬ 
firms the connection between the underdamped rheology 
and the extent of kinetic heating. 
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FIG. 3. Flow curves obtained with a Langevin thermostat 
(black) in an underdamped athermal system and (red) in the 
associated (less underdamped) thermostatted system. 


IV. QUALITATIVE EXPLANATION FOR THE 
NONMONOTONIC FLOW CURVES IN THE 
HIGHLY UNDERDAMPED REGIME 


As can be seen in Fig. 2 of the paper, the flow curves 
at large enough Q are nonmonotonic and thus display a 
minimum (Eim, Sm), where Ei is the rescaled shear rate. 
At this point, using the chain rule. 




9E 


dTK 

OTk 

Eim 

dEi 


as 


Tk=Tk 


(As a technical detail, note that, in practice, the second 
partial derivative - at fixed kinetic temperature - can only 


be calculated in a system with stronger damping ther¬ 
mostatted to Tq = Tk (Eim)). Thus, the presence of a 
minimum results from a competition between two (an¬ 
tagonistic) effects: 

(i) an increase of the sample temperature with the 
shear-rate, which bends the flow curve downwards 

(ii) an intrinsic rate effect at large shear rates, whereby 
the driving gets so strong that particles no longer remain 
at the bottom of the energy wells (collisions, etc.). 

When Q increases, i.e., when the damping weakens, 
the temperature effect (i) gets stronger, so that the mini¬ 
mum is expected to be deeper, which is indeed confirmed 
by Eig. 2 of the main text. 


V. DEVIATIONS FROM THE PREDICTED 

BEHAVIORS OBSERVED AT HIGH KINETIC 
TEMPERATURES 

At high kinetic temperatures TR-(j.e., in particular at 
large shear rates and/or large values of Q), Figs. 2 and 
3 of the main text show that deviations appear from 
the Chattoraj-Lemaitre formula for the temperature- 
dependence of the stress (Eq. 3), and, to a much lower 
extent, from the scaling prediction represented in Fig. 3. 

This is not surprising: Chattoraj and Lemaitre ex¬ 
plicitly derived their for formula in the low-temperature 
limit, and our scaling prediction was based on the hy¬ 
pothesis that regions yield (undergo a plastic event) at 
a local stress Eq, whereas, owing to thermal activation, 
the local stress at yield may be significantly lower at high 
Tk- 
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We study the finite-shear-rate rheology of disordered solids by means of molecular dynamics 
simulations in two dimensions. By systematically varying the damping strength ( in the low- 
temperature limit, we identify two well defined flow regimes, separated by a thin (temperature- 
dependent) crossover region. In the overdamped regime, the athermal rheology is governed by the 
competition between elastic forces and viscous forces, whose ratio gives the Weissenberg number 
Wi oc C 7 ; the macroscopic stress E follows the frequently encountered Herschel-Bulkley law E = 
Eo + ky/wl, with yield stress Eq > 0. In the underdamped (inertial) regime, dramatic changes in 
the rheology are observed for low damping: the flow curve becomes nonmonotonic. This change 
is not caused by longer-lived correlations in the particle dynamics at lower damping; instead, for 
weak dissipation, the sample heats up considerably due to, and in proportion to, the driving. By 
thermostatting more or less underdamped systems, we are able to link quantitatively the rheology 
to the kinetic temperature and the shear rate, rescaled with Einstein’s vibration frequency. 
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Inertia matters in liquid flows. Its presence in the 
Navier-Stokes equations leads to a rich phenomenology 
that vanishes in the overdamped limit of viscous flow. 
However, the effect of damping is rarely heeded (let alone 
analyzed) in the flow of disordered solids, so much so that 
dense colloidal glasses often serve as model systems for 
bulk metallic glasses (BMG) [l|, even though they are 
much more strongly damped. Here, we find that reducing 
the damping can dramatically impact the macroscopic 
rheology. In the inertial regime, the energy input dwells 
longer in the particle momenta before its final dissipa¬ 
tion into the heat bath, thus facilitating plastic flow. We 
provide a quantitative account of this effect in terms of 
simple kinetic heating of the underdamped solid, similar 
to the one observed experimentally during the operation 
of shear bands in BMG ii- 

The damping regime is not the only line of contrast 
among disordered solids: atoms in BMG as well as small 
colloids are heavily influenced by thermal fluctuations 
whereas grains are quasi-athermal; foam bubbles are de¬ 
formable whereas some colloids are close to perfect hard 
spheres. Notwithstanding these contrasting features, vir¬ 
tually all such solids deform similarly, i.e., mostly elas¬ 
tically at small stresses while at larger shear plasticity 
becomes dominant, with a succession of failures of micro¬ 
regions, whose particles rearrange swiftly. These rear¬ 
rangements are triggered W the loading or facilitated 
by thermal activation [ 3 , Q and may interact via the 
long-range elastic deformation that they induce in the 
surrounding medium gill- Based on this generic sce¬ 
nario, multiple simplified rheological models have been 
proposed, generally focusing on the overdamped regime 
|7 hI 2| (nevertheless, the mesoscale elastic response has 


been studied across the damping regimes [Hill). To 
what extent does the presence of inertia alter the pic¬ 
ture? 

In the quasistatic limit, i.e., at vanishing shear rates 7 , 
recent numerical work by Salerno and Robbins has ascer¬ 
tained that the statistics of avalanches fall into distinct 
universality classes in the overdamped vs. underdamped 
regimes (la . 0 - The difference is best illustrated by 
considering the complex, rugged Potential Energy Land¬ 
scape (PEL) in which the system evolves: it climbs up en¬ 
ergy barriers in the phases of elastic loading and abruptly 
slides downhills once the barrier is overcome. For over¬ 
damped systems, this descent suffices to dissipate the en¬ 
ergy stored during loading, while at lower damping the 
inertial force may carry the system over several succes¬ 
sive barriers. This process is then highly directional in 
the PEL and strongly correlated in space and time, which 
renders its modelling quite complex a priori. Some ad hoc 
rules to include inertia in lattice-based models have been 
put forward, such as lowering barriers or yield stresses 
for a certain time after failure (see 0 and references 
in 0 ) and their impact has been emphasized, but the 
validity of these descriptions stands on shaky ground. 

In this Letter, we focus on the steady-state shear 
flow of two-dimensional disordered solids at finite driv¬ 
ing rates and investigate the role of inertia in the vanish¬ 
ing and low temperature limits, with molecular dynamics 
(MD) simulations. We apply simple shear to the binary 
Lennard-Jones glass used in Ref. [1^; it comprises 32,500 
large (type A) particles and 17,500 small (type B) parti¬ 
cles, all of mass m, and has reduced density p = 1.2. The 
equations of motion are based on the Dissipative Particle 
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Dynamics (DPD) scheme and read 




dVjrij) 

9rii 


fi + f, 


D 


( 1 ) 


Here, 


= Ti — r. 


and V{rij) is the interaction po¬ 


tential between particles i and j. The DPD forces fi^’^ 
involve a cutoff function w(r) = 1 —^ if r<rc = 2.5a aa 
and 0 otherwise; = s ^ i^ij) ^ij ^ is a stochas¬ 
tic force, based on the Gaussian white noise Oij (l 8 | 
and due to the coupling to a heat bath maintained at 




temperature Tq and fi^ = y zj j—p—> ij 

is a damping force depending on the relative velocities 
Vij = Vi — Vj. The strength s = 2(kBTo of the cou¬ 
pling to the reservoir depends on the damping strength 
Q and To, and is maintained even if the system departs 
from thermal equilibrium. In the following, C, m, and Tq 
shall be varied, while the particle interactions are kept 
constant. 

The equations of motion, Eqs. [T1 are integrated on 
GPU with the velocity Verlet algorithm. They involve 
forces deriving from four types of stresses: 

(i) the elastic stress, of order S 4 = -^ 4 ^ = 1 , 

• • ■ • n • 

(ii) the viscous stress, of order 777 , where rj fa |l4l | is 
the microscopic viscosity, 

(hi) the inertial pressure, which, in a Bagnold-like pic¬ 
ture [l9j, involves momentum transfers of order maAA'j 
at a rate oc 7 , and is thus proportional to 7717 ^, and 

(iv) the thermal pressure resulting from stochastic 
forces of magnitude ^/C^o- 

Their relative magnitudes are quantified by dimension¬ 
less numbers that characterize the flow regime. In partic¬ 
ular, the importance of viscosity with respect to elasticity 
is measured by the Weissenberg number. 


Wi = Tdissb with Tdiss = 


Sa’ 


and the ratio of inertial over elastic stresses is given by 
Ei^, where 


Ei = Tvibb with Tvib = 

In conjunction with Tq, Wi and Ei fully characterize 
the flow. Nevertheless, to describe the damping regime of 
flow curves, irrespective of the shear rate, it is convenient 
to also introduce 

Ei y'mY.A Tdamp m 

Q = 7777 = -7 - = - with Tdamp = -r; 

Wl C "Tvib C 

if Eq. [1] is assimilated to a damped second-order har¬ 
monic oscillator, Q is the (inertial) quality factor, i.e., 
the number of inertial oscillations in the damping time. 

Onr numerical data confirm the relevance of such di¬ 
mensional analysis: Figures[T]and[5]prove that, at Tg = 0, 




FIG. 1. Athermal flow curves E(Wi,To = 0) in the over¬ 
damped regime Q ;$ 1, for various combinations [C, 7 n,]: 
[C = 1, m = 1] (0), [10,1] (□), [10, 0.1] (A), and [1, 0.1] (V). 
The solid line represents Eq. A flow curve at To = 0.2, 
Q = 1 [1,1] (♦) is also shown. The thin dashed line is a best 
fit to Eq.[3| Inset: E vs. 7 . 


the dependences of the macroscopic shear stress E on C, 
m, and 7 can be condensed into a dependence on the pair 
(Q,Wi), or equivalently (bnt more conveniently when 
Q> 1) (Q,Ei). 

Overdamped dynamics. Let us start by investigating 
the fully overdamped (Brownian or athermal) limit Q —>■ 
0. In the absence of inertia, 7 is best rescaled as Wi. 
At To = 0, the flow curve, plotted in Fig. [U is very well 
described by the Herschel-Bulkley law 

E(Wi,To = 0) = 0.72-h2 Vwi. ( 2 ) 


Interestingly, this description remains very good at fi¬ 
nite values of Q, up to Q Ri 1. Thus, for all Q ^ 1, 
the macroscopic rheology is exclusively governed by the 
competition between elastic and viscous forces. 

Leaving the athermal regime, we observe that imposing 
a finite bath temperature Tg > 0 leads to a decrease of S 
at all shear rates. Regardless of the damping regime, this 
thermal effect is explained by the premature occurrence 
of plastic rearrangements owing to thermal activation: 
the system at Tg > 0 hops over saddle-node points in 
the PEL before the effective potential barriers have com¬ 
pletely flattened under the influence of the driving, which 
interrupts the elastic accumulation of strain; hence the 
lower macroscopic stress [ 2 ^ 


m- 


Perhaps less expectedly, we also find a narrowing of 
the overdamped regime with Tg, that is, the quality fac¬ 
tor Qc(Tg) marking the departure from the scaling with 
Wi decreases with Tg (our data suggest Qc (Tg = 0) ~ 1 
whereas Qc (Tg = 0.2) < 1 but do not allow for greater 
accuracy). A rather general explanation consists in allud¬ 
ing to the excitation of higher-frequency modes at higher 
temperature, these modes having larger specific quality 
factors Q, or to the faster thermalization of the system 
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FIG. 2. Flow curves E(Q,Ei, 0) of athermal underdamped 
systems and E (Q', Ei, Tk (Q, Ei)) of their thermostatted 
counterparts (see text). Symbols are listed in Table |T] Thin 
dashed lines are the best fit to Eq. [31 E = 0.69+ 2 VeI - 

0.17T^/" In ( , where Tk = 0.15Q • Ei. 


(see below). 

Inertial dynamics. On increasing Q, past a small 
crossover region around Qc(To), one enters the under¬ 
damped regime, where the rheology is a priori described 
by the triplet (Q,Ei, Tp). What role does the inertial 
quality factor Q play in that regime? 

In fact, at low damping, Q can no longer be inter¬ 
preted as the number of not-too-damped inertial oscil¬ 
lations within a particle’s cage. Indeed, localized exci¬ 
tations spread in the glass and, owing to nonlinearities, 
thermalize: their energy is redistributed across the whole 
vibrational spectrum. This process occurs over a time Tth 
and expedites the decorrelation of the excitations when 
the velocity damping time '"/c becomes longer than Tth- 
As a result, the velocity autocorrelation functions, which 
reflect single-particle dynamics, gain independence from 
Q, in the quiescent system at Tp > 0 (in Fig. 1 of the 
Supplementary Material (SM), we observe Tth ~ O.lrvib 
at Tp = 0.16). Thus, one is lured into thinking that the 
underdamped rheology is insensitive to Q, in the same 
way as the equilibrium properties of liquids computed 
with MD are independent of the (weak) damping |l8l.l22|. 

Contrary to this thought, the underdamped flow 
curves, plotted in Fig. exhibit dramatic changes at 
large Q (and low Tp), as they become nonmonotonic! 

Clearly, the insensitivity to Q was a fallacy. In fact, 
this parameter also controls energy dissipation in the 
system. When the damping is too weak compared to 
the energy input, the system heats up and strongly de¬ 
parts from thermal equilibrium with the heat reservoir at 
Tp. This is not a numerical artifact: in experiments on 
sheared granular matter, the “granular temperature” dif¬ 
fers from room temperature [23j; temperature rises have 
also been borne out experimentally in shear bands in 


TABLE I. Parameters and symbols used in Figs. [5] and |31 
Inertial quality factors are denoted Q', instead of Q, for the 
thermostatted systems (see text). 



FIG. 3. Shear contribution — Tk — Tq to the kinetic 
temperature measured in underdamped samples vs. Q Ei. 

(V) Data at To = 0.2. The line represents Tj — 0.15 • Q • Ei. 


compressed BMC ii (incidentally, note that a negative 
rate-dependence of the stress, known as “rate-weakening”, 
has also been reported in these materials 0). For “dry” 
systems, heat is actually removed faster in simulations 
than in experiments, where its extraction must proceed 
through the boundaries [l^. Besides, nonmonotonic flow 
curves are not a marginal effect of the DPD thermostat; 
they were also observed by Salerno with a weak Langevin 
thermostat (see Fig. 2.3 of Ref. (^l. 

Taking into account the heating of the sample, we pro¬ 
pose to substitute, in the triplet (Q, Ei, Tp), the reservoir 
temperature Tp with the actual kinetic temperature of 
the sample, Tk = ■ To assess the contribu¬ 

tion Tj of the driving to Tk, we assume that the kinetic 
energy is mainly generated by plastic rearrangements, 
during which the elastic energy ^Sp+y per particle is 
first converted into kinetic energy and then gradually 
dissipated, over a timescale ruamp = ™/c in the under¬ 
damped regime. Thus, the density of simultaneous events 
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is '"/c • "^6 arrive at 


T, 


7 ~ I 


P . Q . Ei. 

Qly 2p 


The scaling law with Q • Ei is in very good agreement 
with the numerical data at Tq = 0, as shown in Fig. [31 
as long as Q ^ 1 (see Section V of SM for a tentative 
rationalization of the slight deviations). The predicted 
(0.29) and measured (0.15) prefactors differ by a factor 
of 2 precisely, probably because the released energy is 
actually equipartitioned between the kinetic and elastic 
degrees of freedom, as for a harmonic oscillator. Further¬ 
more, we observe a Boltzmann distribution, parametrized 
by Ty, of kinetic energies among the particles (Fig. 2 of 
SM), which confirms the status of Ty as the sample tem¬ 
perature. This is consistent with the “quasi-equilibrium” 
situation (at Tk) reported by Xu et al. in strongly- 
sheared athermal systems [^. At finite Tq, we expect 
Tk ~ 2 o -|- Ty, which is entirely compatible with our 
(limited) dataset (Fig. [3]). 

Coming back to the underdamped flow curves. Fig. [2] 
gives ample evidence that the athermal flow curves at 
any Q ^ Qc(0) can be quantitatively reproduced by 
thermostatting a less underdamped (but still inertial) 
system, at Q' < Q, to the shear-rate-dependent tem¬ 
perature TR-(Q,Ei) of the original system; this holds 
true at To > 0 [data at Tq = 0.2 not shown). Put 
differently, E(Q,Ei, Tq) collapses onto a master curve 
S (Ei, Tk (Q, Ei)), irrespective of the value of Q. Thus, Q 
does not impact the underdamped rheology as the inertial 
quality factor, but only via its control of Tk (Q,Ei). It 
follows that inertial vibrations and thermal fluctuations 
have an analogous effect on the rheology: both are “ag¬ 
itation” forces that precipitate rearrangements, but the 
former increase with the shear rate, hence the severe rate¬ 
weakening observed in strongly underdamped systems. 
We should mention that rate-weakening is generally as¬ 
sociated with a flow instability leading to shear-banding 
Hill, but here we have not seen any banding of the 
velocity profiles. We believe that this is due to the ra¬ 
pidity of equilibration through thermal diffusion in small 
systems, which impedes the coexistence of bands sheared 
at different rates, thus (here) at distinct temperatures. 

Chattoraj et al. 12 ill . building on previous work by 
Johnson and Samwer [20|, propounded the following for¬ 
mula for the temperature dependence of the stress. 


/ DyiVe \ 

E(7,To) = E(7,ro = 0)-Ar=/Mn( —, 


(3) 


where A and B are adjustable parameters. Substituting 
To with Tk = Tq + 0.15Q • Ei and 7 with Ei in Eq. [31 we 
obtain predictions in broad agreement with our data, as 
shown in Fig.|3J as long as the flow remains underdamped 
and at low enough Tk- 


These results do not imply that in underdamped sys¬ 
tems inertia can be discarded in favor of temperature. 
Indeed, the collapse onto E (Ei, Tk (Q, Ei)) breaks down 
for Q < Qc(ro), which highlights the operativeness of an 
inertial mechanism at Q > Qo(To)) responsible e.g. for 
the scaling of the inverse attempt frequency (multiplied 
by 7 ) with Ei, and not Wi. Still, it is noteworthy that the 
collapse holds down to values of Q in the crossover region; 
in particular, systems at Q = 1 display a macroscopic 
rheology close to the fully overdamped one at Tq = 0, 
while a scan through their higher-temperature response 
gives access to the strongly underdamped rheology. 

In summary, the variations of the macroscopic rheology 
of a model disordered solid with damping strength () (or 
particle mass m) can be collapsed into two flow regimes. 
When Q = Ei/Wi is smaller than a threshold Qc(7o), the 
system is overdamped. It is widely accepted that foams, 
concentrated emulsions, and dense colloidal suspensions 
belong in this regime. At fixed Tq, in particular Tq = 0, 
the flow curves only depend on Wi, which proves that the 
competition between the elastic interactions imposed by 
the PEL and dissipation forces dominates the rheology 
of these systems. This is compatible with the rheological 
models proposed by us and others in [iSlIIlMlH , but 


rules out all explanations based on the transverse sound 
velocity Cg (which affects Ei and Q, but not Wi). 

Such explanatory scenarios based on Cg could be 
valid in the (moderately) underdamped regime, at Q > 
Qr(Tn). As a noteworthy example, Lemaitre and Carol! 
[33 suggested the following scenario, later taken up and 
revised in [H, [s^l : avalanches of plastic events spread 
at speed Cg and their spreading is limited by the driving, 
which generates independent plastic events. The ensuing 
incomplete plastic relaxation explains the increasing flow 
curve. However, the athermal MD system used in 
appears similar to ours with Q Ri 0.2 < Qc(Tb — 0)- 

For even more strongly underdamped systems, at Q ^ 
10^, the flow curve becomes nonmonotonic at low bath 
temperature. Surprisingly, this transition has never been 
analyzed before, although the threshold for Q does not 
seem unrealistically large: a crude estimate for a suspen¬ 
sion of frictionless grains (of density p and radius a) in 
a solvent of viscosity 77 gives Q r; o ^ showed 

that variations in the inertial properties of the material 
played no role per se in the transition; instead, the lat¬ 
ter originates from the insufficient energy dissipation at 
large Q, which causes the sample to heat up (and hence, 
relax stress) all the more as the driving is fast, with the 
scaling law Ty oc Q • Ei. 

This rate-weakening mechanism persists until it is 
counterbalanced by the standard collisional increase of E 
at high rates; this results in a minimum in the flow curve 
(see Section IV of SM for a discussion). The mechanism 
is reminiscent of the one producing a shear-banding insta¬ 
bility in the Soft Glassy Rheology (SGR)’s variant pro¬ 
posed by Fielding et al. [^. In SGR, material subunits 






5 


possess a (widely distributed) energy barrier for yield¬ 
ing, which decreases as the material is loaded. Yielding 
is then activated by an effective mechanical temperature 
X. In the variant of Ref. (s^, x is coupled to the local 
plastic activity and thus increases with the shear rate. 
In a similar fashion, in the Shear Transformation Zone 
theory, the strain may localize via a coupling between 
the strain rate and the “configurational disorder temper¬ 
ature” [ 3 ^. The major conceptual divergence between 
these approaches and our observations in severely under¬ 
damped systems is the (effective or kinetic) nature of the 
temperature. 

This difference echoes a vast debate in the metallic 
glass community regarding the origin of the softening of 
shear bands: does the band persist by softening because 
of heat production, hence, higher local temperatures, or, 
perhaps more probably, because of local configurational 
changes (in free volume or density), while the tempera¬ 
ture rise is but a side-effect ii [ 25 ! ? Our findings do 
not contribute to settling this question, but they do cer¬ 
tainly call for a clarification of the description of damping 
in rheological models. 
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